
# generate appendix C3, figures C1 to C8

library(rio) # to load data
library(tidyverse) # to reshape data
library(forecast) # to run time series diagnostics

# set your working directory
setwd("~/replication_files")

monthly_data <- import("data/monthly_data.csv")

# to test for stationarity, need to use complete observations
ts_data <- monthly_data %>%
  select(iso3c,time,log_sovereign_issued_1yr_amount,log_oil_production,terms_of_trade,resource_rents_perc) %>%
  drop_na()

# confidence intervals
df_ci <- ts_data %>% 
group_by(iso3c) %>% 
  summarise(ci = qnorm((1 + 0.95)/2)/sqrt(n()))

# figure C1
df_acf1 <- ts_data %>% 
  group_by(iso3c) %>% 
  summarise(list_acf = list(acf(log_sovereign_issued_1yr_amount, plot = FALSE))) %>%
  mutate(acf_vals = purrr::map(list_acf, ~as.numeric(.x$acf))) %>% 
  select(-list_acf) %>% 
  unnest(cols = c(acf_vals)) %>% 
  group_by(iso3c) %>% 
  mutate(lag = row_number() - 1)

ggplot(df_acf1, aes(x = lag, y = acf_vals)) +
  geom_bar(stat = "identity", width=.3) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color = "blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color = "blue", linetype="dotted") +
  labs(x="Lag", y ="ACF") +
  facet_wrap(~iso3c, ncol = 5) + theme_classic() 

ggsave("figures/figureC1.pdf", height = 10, width = 8)

# figure C2
df_pacf1 <- ts_data %>% 
  group_by(iso3c) %>% 
  summarise(list_pacf = list(pacf(log_sovereign_issued_1yr_amount, plot = FALSE))) %>%
  mutate(pacf_vals = purrr::map(list_pacf, ~as.numeric(.x$acf))) %>% 
  select(-list_pacf) %>% 
  unnest(cols = c(pacf_vals)) %>% 
  group_by(iso3c) %>% 
  mutate(lag = row_number() - 1)

ggplot(df_pacf1, aes(x = lag, y = pacf_vals)) +
  geom_bar(stat = "identity", width=.3) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color = "blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color = "blue", linetype="dotted") +
  labs(x="Lag", y ="PACF") +
  facet_wrap(~iso3c, ncol = 5) + theme_classic() 

ggsave("figures/figureC2.pdf", height = 10, width = 8)

# figure C3
df_acf2 <- ts_data %>% 
  group_by(iso3c) %>% 
  summarise(list_acf = list(acf(resource_rents_perc, plot = FALSE))) %>%
  mutate(acf_vals = purrr::map(list_acf, ~as.numeric(.x$acf))) %>% 
  select(-list_acf) %>% 
  unnest(cols = c(acf_vals)) %>% 
  group_by(iso3c) %>% 
  mutate(lag = row_number() - 1)

ggplot(df_acf2, aes(x = lag, y = acf_vals)) +
  geom_bar(stat = "identity", width=.3) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color = "blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color = "blue", linetype="dotted") +
  labs(x="Lag", y ="ACF") +
  facet_wrap(~iso3c, ncol = 5) + theme_classic() 

ggsave("figures/figureC3.pdf", height = 10, width = 8)

# figure C4
df_pacf2 <- ts_data %>% 
  group_by(iso3c) %>% 
  summarise(list_pacf = list(pacf(resource_rents_perc, plot = FALSE))) %>%
  mutate(pacf_vals = purrr::map(list_pacf, ~as.numeric(.x$acf))) %>% 
  select(-list_pacf) %>% 
  unnest(cols = c(pacf_vals)) %>% 
  group_by(iso3c) %>% 
  mutate(lag = row_number() - 1)

ggplot(df_pacf2, aes(x = lag, y = pacf_vals)) +
  geom_bar(stat = "identity", width=.3) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color = "blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color = "blue", linetype="dotted") +
  labs(x="Lag", y ="PACF") +
  facet_wrap(~iso3c, ncol = 5) + theme_classic() 

ggsave("figures/figureC4.pdf", height = 10, width = 8)

# figure C5
df_acf3 <- ts_data %>% 
  group_by(iso3c) %>% 
  summarise(list_acf = list(acf(log_oil_production, plot = FALSE))) %>%
  mutate(acf_vals = purrr::map(list_acf, ~as.numeric(.x$acf))) %>% 
  select(-list_acf) %>% 
  unnest(cols = c(acf_vals)) %>% 
  group_by(iso3c) %>% 
  mutate(lag = row_number() - 1)

ggplot(df_acf3, aes(x = lag, y = acf_vals)) +
  geom_bar(stat = "identity", width=.3) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color = "blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color = "blue", linetype="dotted") +
  labs(x="Lag", y ="ACF") +
  facet_wrap(~iso3c, ncol = 5) + theme_classic() 

ggsave("figures/figureC5.pdf", height = 10, width = 8)

# figure C6
df_pacf3 <- ts_data %>% 
  group_by(iso3c) %>% 
  summarise(list_pacf = list(pacf(log_oil_production, plot = FALSE))) %>%
  mutate(pacf_vals = purrr::map(list_pacf, ~as.numeric(.x$acf))) %>% 
  select(-list_pacf) %>% 
  unnest(cols = c(pacf_vals)) %>% 
  group_by(iso3c) %>% 
  mutate(lag = row_number() - 1)

ggplot(df_pacf3, aes(x = lag, y = pacf_vals)) +
  geom_bar(stat = "identity", width=.3) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color = "blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color = "blue", linetype="dotted") +
  labs(x="Lag", y ="PACF") +
  facet_wrap(~iso3c, ncol = 5) + theme_classic() 

ggsave("figures/figureC6.pdf", height = 10, width = 8)

# figure C7
df_acf4 <- ts_data %>% 
  group_by(iso3c) %>% 
  summarise(list_acf = list(acf(terms_of_trade, plot = FALSE))) %>%
  mutate(acf_vals = purrr::map(list_acf, ~as.numeric(.x$acf))) %>% 
  select(-list_acf) %>% 
  unnest(cols = c(acf_vals)) %>% 
  group_by(iso3c) %>% 
  mutate(lag = row_number() - 1)

ggplot(df_acf4, aes(x = lag, y = acf_vals)) +
  geom_bar(stat = "identity", width=.3) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color = "blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color = "blue", linetype="dotted") +
  labs(x="Lag", y ="ACF") +
  facet_wrap(~iso3c, ncol = 5) + theme_classic() 

ggsave("figures/figureC7.pdf", height = 10, width = 8)

# figure C8
df_pacf4 <- ts_data %>% 
  group_by(iso3c) %>% 
  summarise(list_pacf = list(pacf(terms_of_trade, plot = FALSE))) %>%
  mutate(pacf_vals = purrr::map(list_pacf, ~as.numeric(.x$acf))) %>% 
  select(-list_pacf) %>% 
  unnest(cols = c(pacf_vals)) %>% 
  group_by(iso3c) %>% 
  mutate(lag = row_number() - 1)

ggplot(df_pacf4, aes(x = lag, y = pacf_vals)) +
  geom_bar(stat = "identity", width=.3) +
  geom_hline(yintercept = 0) +
  geom_hline(data = df_ci, aes(yintercept = -ci), color = "blue", linetype="dotted") +
  geom_hline(data = df_ci, aes(yintercept = ci), color = "blue", linetype="dotted") +
  labs(x="Lag", y ="PACF") +
  facet_wrap(~iso3c, ncol = 5) + theme_classic() 

ggsave("figures/figureC8.pdf", height = 10, width = 8)


## bonus: test for stationarity, by variable and by country

# dv: ln sovereign amount issued
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "ARG"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,1)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "BOL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,2)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "BRA"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,1) with drift
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "CHL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "COL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,2)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "CRI"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "DOM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA (0,1,2)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "ECU"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,0,1) with non-zero mean
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "GTM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,0,0)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "GUY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,0,0)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "HND"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,0,1) with non-zero mean
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "JAM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,4)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "MEX"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "NIC"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,1)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "PAN"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,0,0)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "PER"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "PRY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "SUR"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(5,1,2)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "SLV"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(2,1,3)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "TTO"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,3)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "URY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,3)
auto.arima(ts_data$log_sovereign_issued_1yr_amount[ts_data$iso3c == "VEN"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)

# iv: resource rents
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "ARG"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "BOL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "BRA"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "COL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "CRI"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "DOM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "ECU"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "GTM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "GUY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "HND"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "JAM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "MEX"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "PAN"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "PER"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "SLV"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$resource_rents_perc[ts_data$iso3c == "VEN"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)

# iv: ln oil and gas production
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "ARG"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(3,1,4)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "BOL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,3)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "BRA"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "CHL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "COL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "CRI"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "DOM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # no production
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "ECU"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,2)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "GTM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,3)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "GUY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "HND"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "JAM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "MEX"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(2,2,2)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "NIC"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "PAN"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "PER"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,2)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "PRY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "SUR"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "SLV"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "TTO"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,2) with drift
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "URY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,0)
auto.arima(ts_data$log_oil_production[ts_data$iso3c == "VEN"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(4,1,0)


# iv: commodity price index
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "ARG"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,0)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "BOL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,1)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "BRA"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,3)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "CHL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,0)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "COL"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,2)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "CRI"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(4,1,2)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "DOM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "ECU"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,3)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "GTM"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(3,1,1)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "GUY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(3,1,0)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "HND"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(4,1,2)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "JAM"], 
           test=c("kpss","adf","pp"), trace=TRUE)  # ARIMA(2,1,1)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "MEX"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(2,1,1)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "NIC"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,2)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "PAN"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(2,1,4)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "PER"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,0)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "PRY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(0,1,1)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "SLV"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(2,1,1)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "SUR"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(4,1,1 with drift)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "TTO"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(2,1,0)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "URY"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(1,1,4)
auto.arima(ts_data$terms_of_trade[ts_data$iso3c == "VEN"], 
           test=c("kpss","adf","pp"), trace=TRUE) # ARIMA(3,1,1)
